﻿using System;
using System.Collections.Generic;
using Allegro.Integrator;

namespace Allegro.Integrator.RungeKutta
{
    /// <summary>
    /// This is a stage 8, order (5,6) Runge-Kutta pair published by Jim Verner.
    /// Documentation of this (5,6) Runge-Kutta pair:
    /// http://www.math.sfu.ca/~jverner/New56Pair.040729.pdf
    /// </summary>
    public class RKTableau_56_Verner02 : ButcherTableau
    {
        public RKTableau_56_Verner02() : base(8)
        {
            // Initialize the vectors
            c[0] = 0.0;
            c[1] = 1.0 / 7.0;
            c[2] = 2.0/9.0;
            c[3] = 3.0/7.0;
            c[4] = 2.0/3.0;
            c[5] = 3.0/4.0;
            c[6] = 1.0;
            c[7] = 1.0;

            a[1, 0] = 1.0/7.0;

            a[2, 0] = 4.0 / 81.0;
            a[2, 1] = 14.0 / 81.0;

            a[3, 0] = 291.0 / 1372.0;
            a[3, 1] = -27.0 / 49.0;
            a[3, 2] = 1053.0 / 1372.0;

            a[4, 0] = 86.0 / 297.0;
            a[4, 1] = -14.0 / 33.0;
            a[4, 2] = 42.0 / 143.0;
            a[4, 3] = 1960.0 / 3861.0;

            a[5, 0] = -267.0 / 22528.0;
            a[5, 1] = 189.0 / 704.0;
            a[5, 2] = 63099.0 / 585728.0;
            a[5, 3] = 58653.0 / 366080.0;
            a[5, 4] = 4617.0 / 20480.0;

            a[6, 0] = 10949.0 / 6912.0;
            a[6, 1] = -69.0 / 32.0;
            a[6, 2] = -90891.0 / 68096.0;
            a[6, 3] = 112931.0 / 25920.0;
            a[6, 4] = -69861.0 / 17920.0;
            a[6, 5] = 26378.0 / 10773.0;

            a[7, 0] = 1501.0 / 19008.0;
            a[7, 1] = -21.0 / 88.0;
            a[7, 2] = 219519.0 / 347776.0;
            a[7, 3] = 163807.0 / 926640.0;
            a[7, 4] = -417.0 / 640.0;
            a[7, 5] = 1544.0 / 1539.0;
            a[7, 6] = 0.0;

            // High order weights
            b[0] = 79.0 / 1080.0;
            b[1] = 0.0;
            b[2] = 19683.0 / 69160.0;
            b[3] = 16807.0 / 84240.0;
            b[4] = 0.0;
            b[5] = 2816.0 / 7695.0;
            b[6] = 0.01;
            b[7] = 187.0 / 2800.0;

            // Low order weights
            bh[0] = 763.0 / 10800.0;
            bh[1] = 0.0;
            bh[2] = 59049.0 / 197600.0;
            bh[3] = 88837.0 / 526500.0;
            bh[4] = 243.0 / 4000.0;
            bh[5] = 12352.0 / 38475.0;
            bh[6] = 0.0;
            bh[7] = 2.0 / 25.0;
        }
    }
}
